How to create a marker map with leaflet

Resources (to work with spatial vector data in R)

You need R, RStudio, and some packages to get started. If you are completely new, I recommend resources from RLadies Sydney (BasicBasics and CleanItUp units are must-read).

In R, Working with tabular data using dplyr package is straightforward, |>or%>% pipe operators and these functions: mutate, select, filter and group_by along with 4 types of join should get you far.

The most popular package to work with spatial data in R is sf. An sf object extends regular tabular data frames to include a geometry column, which contains the spatial information. Thus, a simple explanation on how to read shapefiles with sf should be enough to get you started. Then, those mentioned dplyr functions will come in handy.

Motivation

To build maps in R, the two most popular packages are mapview and leaflet. For this report, I chose leaflet due to its flexibility and customizability.

Here’s the reasoning behind my choice and map design:

• To show transportation hubs and major roads → Map tiles
• To display warehouse locations with details → Markers with popups
• To distinguish warehouse sizes → Marker size
• To indicate warehouse type (3PL vs. in-house) → Marker color
• To compare pre- and post-COVID landscapes → Two separate layers: before and after
• To provide population context → Sub-layer of population density dots (since this is not the main focus)

Execution

1. Preparing data

Fullfillment Centers data

Data on fulfillment centers in Greater Melbourne was collected through manual searches on Google Maps, company website visits, news articles, and LinkedIn posts for 8 piece of information:

• Name;

• Address;

• Since (the year the warehouse was established).

• Type (in-house or 3PL);

• Area (size of the warehouse, categorized as Small, Medium, Large);

• Note (any additional information about the warehouse).

• lat: (latitude coordinate of address);

• lon (longitude coordinate of address);

lon and lat were obtained using the geocode function which retrieves latitude and longitude based on the address. Data is then turned into a sf object using sf::st_as_sf(coords = c("lon", "lat")

import requests

def geocode(address):
    url = "https://nominatim.openstreetmap.org/search"
    params = {
        'q': address,
        'format': 'json'
    }
    response = requests.get(url, params=params, headers={'User-Agent': 'Your App Name'})
    if response.status_code == 200:
        data = response.json()
        if data:
            return data[0]['lat'], data[0]['lon']
    return None, None

Data is available for download from fullfillmentCenter_sf.zip

fullfillmentCenter_sf

Simple feature collection with 28 features and 6 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 144.729 ymin: -38.05242 xmax: 145.2347 ymax: -37.578
Geodetic CRS:  GCS_GDA2020
# A tibble: 28 × 7
   Name                Address Type  Since Area  Note              geometry
   <chr>               <chr>   <chr> <dbl> <chr> <chr>          <POINT [°]>
 1 Officeworks Distri… 75 Hor… inho…  2023 Large Tota… (144.8108 -37.84454)
 2 GPC Melbourne Dist… 24 Wil… inho…  2023 Large Tota…  (144.948 -37.67748)
 3 Aldi - Distributio… 60 Swa… inho…  2010 Mega  Tota… (144.7756 -37.80179)
 4 Myer National Dist… 73 Mom… inho…  2023 Large Tota…  (144.729 -37.76164)
 5 Somerton Distribut… 20-50 … 3PL    2010 Mega  7134… (144.9574 -37.64112)
 6 Kmart Distribution… 2-12 B… inho…  2018 Mega  Tota… (144.7609 -37.84188)
 7 Coles Distribution… 484 Do… inho…  2024 Mega  Cole… (144.7338 -37.82371)
 8 Target Distributio… 15 Nat… inho…  2017 Mega  Targ… (144.7589 -37.83151)
 9 John Deere Limited… 1 John… inho…  2020 Large Tota… (144.8059 -37.81015)
10 JB Hi-Fi Melbourne… 30 Aga… inho…  2020 Large 3048… (144.7464 -37.79908)
# ℹ 18 more rows

Population Density data

How to generate dots density data is explained in the next section.

Data is available for download from points_sf.zip

points_sf

Simple feature collection with 57470 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 141.0321 ymin: -38.87474 xmax: 149.7605 ymax: -34.12203
Geodetic CRS:  GCS_GDA2020
# A tibble: 57,470 × 2
     FID             geometry
   <dbl>          <POINT [°]>
 1     0  (143.802 -37.55391)
 2     1  (143.804 -37.55111)
 3     2 (143.8023 -37.54956)
 4     3 (143.8039 -37.55215)
 5     4 (143.7962 -37.55324)
 6     5 (143.8123 -37.56792)
 7     6 (143.8129 -37.56038)
 8     7 (143.8142 -37.56515)
 9     8 (143.8009 -37.55428)
10     9 (143.7993 -37.55614)
# ℹ 57,460 more rows

2. Recreating the map

map <- leaflet::leaflet() |> 
  
1  leaflet::addProviderTiles("OpenStreetMap") |>
  
2  leaflet::fitBounds(
    lng1 = 144.72903, lat1 = -38.05242,
    lng2 = 145.23470, lat2 = -37.57800
  ) |>

7  leafgl::addGlPoints(
    data = points_sf,
    fillColor = "orange",
    fillOpacity = 0.5,
    radius = 2,
    group = "Population Density Dots"
  ) |>
  
3  leaflet::addCircleMarkers(
    data = fullfillmentCenter_sf,
    group = "After COVID",
    radius = ~dplyr::case_when(
      Area == "Small" ~ 4,
      Area == "Medium" ~ 6,
      Area == "Large" ~ 8,
      TRUE ~ 10
    ),
    color = ~dplyr::case_when(
      Type == "inhouse" ~ "red",
      Type == "3PL" ~ "navy",
      TRUE ~ "grey"
    ),
    stroke = TRUE,
    weight = 0.1,
    fillOpacity = 0.5,
    label = ~Name,
    popup = ~paste0(
      "<b>Name:</b> ", Name, "<br>",
      "<b>Size:</b> ", Area, "<br>",
      "<b>Since:</b> ", Since, "<br>",
      "<b>Note:</b> ", Note, "<br>"
    )
  ) |> 
  
  leaflet::addCircleMarkers(
    data = fullfillmentCenter_sf %>% dplyr::filter(Since < 2020),
    group = "Before COVID",
    radius = ~dplyr::case_when(
      Area == "Small" ~ 4,
      Area == "Medium" ~ 6,
      Area == "Large" ~ 8,
      TRUE ~ 10
    ),
    color = ~dplyr::case_when(
      Type == "inhouse" ~ "red",
      Type == "3PL" ~ "navy",
      TRUE ~ "grey"
    ),
    stroke = TRUE,
    weight = 0.1,
    fillOpacity = 0.5,
    label = ~Name,
    popup = ~paste0(
      "<b>Name:</b> ", Name, "<br>",
      "<b>Size:</b> ", Area, "<br>",
      "<b>Since:</b> ", Since, "<br>",
      "<b>Note:</b> ", Note, "<br>"
    )
  ) |> 
  
4  leaflet::addLegend(
    position = "bottomright",
    colors = c("Inhouse" = "red", "3PL" = "navy"),
    labels = c("Inhouse", "3PL"),
    title = "Warehouse Type"
  ) |>
  
5  leaflet::addControl(
    html = "<div style='background:white; padding:4px; font-size:12px; border-radius:4px; display:none;'>
              <span style='display:inline-block; width:10px; height:10px; background:orange; border-radius:50%; margin-right:4px;'></span> = 100 persons
            </div>",
    position = "bottomleft"
  ) |>
  
6  leaflet::addLayersControl(
    baseGroups = c("After COVID", "Before COVID"),
    overlayGroups = c("Population Density Dots"),
    options = layersControlOptions(collapsed = FALSE)
  )

1

Base Map Tiles
"OpenStreetMap" was used as the base tile. You can explore more tile options by changing the provider name in addProviderTiles().

2

Initial Map View
fitBounds() sets the initial view using the bounding box of Greater Melbourne.

3

Warehouse Markers
addCircleMarkers() adds two sets of warehouse markers (before and after COVID), customized by: popup for details on click, color for warehouse type (e.g., 3PL vs. in-house), radius for warehouse size

4

Legend addLegend() displays a legend for warehouse types.

5

Population Density Key
addControl() adds a toggleable key explaining the population density dots. It appears only when the corresponding layer is active.

6

Layer Control Panel
addLayersControl() enables toggling between: Pre-COVID warehouses, Post-COVID warehouses, Population density dots

7

Population Density Dots
addGlPoints() from the leafgl package is used to efficiently render thousands of population dots, providing better performance than addCircleMarkers() for large datasets.